Classical simulation of the Hubbard-Holstein dynamics with optical waveguide lattices 
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A classical analog simulator of the two-site Hubbard-Holstein model, describing the dynamics 
of two correlated electrons coupled with local phonons, is proposed based on light transport in 
engineered optical waveguide arrays. Our photonic analog simulator enables to map the temporal 
dynamics of the quantum system in Fock space into spatial propagation of classical light waves in 
the evanescently-coupled waveguides of the array. In particular, in the strong correlation regime the 
periodic temporal dynamics, related to the excitation of Holstein polarons with equal energy spacing, 
can be visualized as a periodic self-imaging phenomenon of the light beam along the waveguide 
array and explained in terms of generalized Bloch oscillations of a single particle in a semi-infinite 
inhomogeneous tight-binding lattice. 
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I. INTRODUCTION 



In recent years, there has been a renewed and increas- 
ing interest in the development of quantum analog sim- 
ulators, i.e. controllable quantum systems that can be 
used to imitate other quantum systems (see, for instance, 
[l| and references therein). Quantum analog simulators 
would be very useful for a wide variety of problems in 
physics, chemistry, and biology. In particular, they have 
proven to be able to tackle complex quantum many- 
body problems in condensed-matter physics ^2ni], such 
as correlated electrons or quantum magnetism. Quan- 
tum analog simulators are also useful to test predic- 
tions or states of matter difficult to access experimen- 
tally in the originally proposed system, such as dynam- 
ical and extreme conditions [1]. Recent realizations of 
quantum analog simulators of many-body problems have 
been based mainly on atoms 0, , trapped ions 0, @- 
[iQil , nuclear magnetic resonance and single photons 
[H, [l3| . As compared to quantum simulators based on 
atoms or ions, those employing single photons offer the 
feasibility to individually address and access the dynam- 
ics of single particles, and are thus suited to simulate 
the physics of few interacting particles [13]. Unlike elec- 
trons, photons do not interact with each other, espe- 
cially at the few photon level, and hence for the sim- 
ulation of quantum many-body problems of condensed- 
matter physics particle interaction has to be introduced 
in a suitable manner. At the few photon level, inter- 
action is generally mimicked using purely linear optical 
systems by measurement-induced nonlinearity [1^. In 
this way, simplified toy models of many-body condensed 
matter physics have been simulated, such as frustrated 
Heisenberg dynamics in a spin- 1/2 tetramer [l3|. On the 
other hand, propagation of classical light in waveguide- 
based optical lattices has provided an experimentally ac- 
cessible test bed to mimic in a purely classical setting 
single-particle coherent phenomena of solid-state physics 
fl5, 16], such as Bloch oscillations p!?!, d ynamic local- 
ization [18], and Anderson localization [19| to mention a 
few. Recently, the possibility to simulate Bose-Hubbard 



and Fermi-Hubbard models of few interacting bosons or 
fermions in coupled-waveguide linear optical structures 
using classical light beams has been proposed as well 
(2Q| , [2H . The main idea is that the temporal evolution 
of a few interacting particles in Fock space can be con- 
veniently mapped into linear spatial propagation of dis- 
cretized light [22| in engineered evanescently-coupled op- 
tical waveguides. As compared to quantum simulators 
using single photons, atoms or trapped ions, such a kind 
of classical simulators of quantum many-body problems 
using classical light would benefit for the much simpler 
experimental implementation, the rather unique possi- 
bility to provide a direct access and visualization of the 
temporal dynamics of the quantum system in Fock space, 
and the ability to prepare the system in a rather arbitrary 
(for example highly entangled) state. 
In condensed-matter, strong correlation effects and local- 
ization can occur in metallic systems due both to strong 
electron-electron(e-e) interactions and strong electron- 
phonon (e-ph) coupling. Evidences of strong e-ph inter- 
actions have been reported in such important materials 
as cuprates, fullerides, and manganites (see, e.g., [23] and 
references therein). The interplay of e-e and e-ph interac- 
tions in these correlated systems leads to coexistence of or 
competition between various phases such as superconduc- 
tivity, charge-density- wave or spin-density- wave phases, 
or formation of novel non- Fermi liquid phases, polarons, 
bipolarons, etc. One of the simplest theoretical model 
that accounts for the interplay between these two types 
of interactions is the Holstein-Hubbard (HH) model [2J- 
i26]. Even though the HH model provides an oversim- 
plified description of both e-e and e-ph interactions, it 
retains what are thought to be the relevant ingredients 
of a system in which electrons experience simultaneously 
an instantaneous short-range repulsion and a phonon- 
mediated retarded attraction. Actually, in spite of its 
formal simplicity, it is not exactly solvable even in one 
dimension. As quantum simulators of Fermi-Hubbard 
models, based on either fermionic atoms in optical lat- 
tices [3I, [23 or electrons in artificial quantum-dot crystals 
[28], have attracted great interest, analog simulators of 
HH models have not received special attention yet. The 
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smallest configuration for studying and simulating the 
HH dynamics is the dimer HH model, which describes 
two correlated electrons hopping between two adjacent 
sites [26]. In In this work, a classical analog simulator is 
theoretically proposed for the two-site HH model, based 
on light transport in engineered waveguide lattices, which 
is capable of reproducing the temporal dynamics of the 
quantum model in Fock space as spatial beam evolution 
in the photonic lattice. The paper is organized as follows. 
In Sec. II, the HH model is briefly reviewed, and the dy- 
namical equations of Fock-state amplitudes are derived. 
In Sec. HI, a classical analog simulator based on light 
transport in engineered waveguide arrays is described, 
and numerical results are presented for waveguide struc- 
tures that simulate the HH model for different electron 
correlation regimes. In particular, in the strong corre- 
lation regime the quasi-periodic temporal dynamics, re- 
lated to the excitation of Holstein polarons, is visualized 
as a periodic self-imaging phenomenon of the light beam 
that propagates along the waveguide array and explained 
in terms of generalized Bloch oscillations of a single par- 
ticle in an inhomogeneous biased tight-binding lattice. 
Finally, in Sec. IV the main conclusions are outlined. 



II. THE HUBBARD-HOLSTEIN MODEL: 
FOCK-SPACE REPRESENTATION 

A. The two-site Hubbard-Holstein model 

The two-site HH model describes two correlated elec- 
trons hopping between two adjacent sites (a diatomic 
molecule), each of which exhibits an optical mode with 
frequency Q. This model has been often used in 
condensed-matter physics as a simplified toy model to 
describe some basic features of polaron and bipolaron dy- 
namics (see, for instance, [26]). The two-site HH Hamil- 
tonian can be separated into two terms [26]. One de- 
scribes a shifted oscillator that does not couple to the 
electronic degrees of freedom, and can be thus disre- 
garded. The other describes the effective e-ph system 
where phonons couple directly with the electronic degrees 
of freedom. The corresponding Hamiltonian, with h = 1, 
reads explicitly [26] 



H = He -\- Hph + Ve-ph 



(1) 



where 



+ U (ni,^ni,; + n2,t^2,i) (2) 
Hph = ^b^b (3) 

(4) 



Ve-ph = (ni - n2) (b + b"^^ . 



V2 

In the previous equations, cJa (^i,cr) are the fermionic 
creation (annihilation) operators for the electron at site 



i with spin a (i = 1,2, a =t,i), ^i,cr = 4cr^*,cr and 
= ^i,t + are the electron occupation numbers, g is 
the on-site e-ph coupling strength, t is the hopping am- 
plitude between adjacent sites, U is the on-site Coulomb 
interaction energy, and (b) is creation (annihilation) 
bosonic operator for the oscillator. For the two-site two- 
electron system there are six electronic states. Three of 
these states are degenerate, with zero energy He\il^Ti) = 
{i = 1, 2, 3), and belong to the triplet states as follows: 

\M = cl,T4,tlO)e (5) 

I^Tl) = ^ (cl.fcL + cU4,t) I0)e (6) 

IVts) = cl44,JO)e. (7) 

The triplet states are not coupled with the 6-oscillator, 
and can be thus disregarded. The three other eigenstates 
of He^ denoted by |— ), \s-\-) and |5— ), are constructed 
from the singlet states and read explicitly 



|0)e 



= sin(9|s) - cos6>|+) 
\s-) = cose\s) +sin6>|+) 



where we have set 
and 




U 



(8) 

(9) 
(10) 
(11) 



(12) 
(13) 

(14) 



For these states one has He\ — ) = U\ — ) and He\s±) 
E±\s±), where 



E± = 



U ± Vf/2 + 16^2 



(15) 



The e-ph interaction term Ve-ph couples the electronic 
eigenstates | — ) and of He with the phonon states 
\n)ph = {l/V^P''\0)ph of Hph. 

B. Fock-space representation 

In order to realize a photonic analog simulator of the 
HH Hamiltonian, it is worth expanding the state vector 
\^p{t)) of the electron-phonon system in Fock space as 
follows 



l^t) = E {fnit)\-) + 9n\s+) - hn{t)\s-)] (g) |n)j,^. 

(16) 



n=0 
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The evolution equations for the amphtude probabihties 
fn^ 9n and hn in Fock space, as obtained from the 
Schrodinger equation idt\il){t)) = H\il){t))^ read exphc- 
itly 

- (17) 
= (^+ + ^^)9n - l^nfn-1 - ^n+l/n+1 (18) 
= (^- + nn)hn - pnfn-1 - Pn+l/n+1 (19) 

(n = 0, 1, 2, 3, ...), where we have set 

tvn = V2ngcos0 , pn = V2ngsm0 (20) 

and where we have assumed /_i = g^i = h-i = for 
definiteness. 



III. PHOTONIC ANALOG SIMULATOR 
A. Waveguide array design 

Equations (17-19) govern the temporal dynamics in 
Fock space of the HH Hamiltonian, and their op- 
tical implementation can be achieved using a sim- 
ple linear optical network because the coupled equa- 
tions of the Fock-state amplitudes are linear equations. 
Let us first notice that Eqs. (17-19) decouple into two 
sets of equations involving the former the amplitudes 
{/2n(^),^2n+i(^),^2n+i(^)}, the latter the amplitudes 
{/2n+i(^),^2n(^),^2n(^)} (n = 0, 1 , 2, 3, . . .) • A Schematic 
of the optical networks that provide analog simulators of 
the two sets of equations are depicted in Fig.l. Their 
physical implementation can be obtained by using either 
arrays of coupled optical waveguides with engineered cou- 
pling rates and propagation constants, where the tempo- 
ral evolution of the Fock-state amplitudes is mapped into 
spatial propagation of the light fields along the axis of the 
various waveguides, or arrays of coupled optical cavities 
with engineered coupling rates and resonance frequencies, 
in which the temporal evolution of the Fock-state ampli- 
tudes is mimicked by the temporal evolution of the light 
field stored in the various cavities. Owing to the current 
technological advances in manufacturing engineered ar- 
rayed waveguides and related imaging techniques based 
on femtosecond laser writing (see, for instance, [16] and 
references therein), we will consider here specifically an 
optical analog simulator of the HH model based on prop- 
agation of monochromatic light waves in coupled optical 
waveguides. Moreover, since the dynamics of the two 
networks of Fig.l is decoupled, without loss of general- 
ity we can limit to consider initial conditions that excite 
only one of them, for example the one corresponding to 
Fig. 1(a). 

Light transport of a classical monochromatic field at 
wavelength A in a weakly-guiding arrays of coupled op- 
tical waveguides with a straight optical axis in the 



geometrical setting schematically depicted in Fig. 1(a), is 
described by the optical Schodinger equation for the elec- 
tric field envelope (j){x^ z) [15] 

where A = A/(27r) is the reduced wavelength of photons. 
Us is the refractive index of the substrate at wavelength 
A, V{x^ y) ns — n{x, y) is the so-called optical potential, 
and n{x^y) is the refractive index profile of the optical 
structure defining the various waveguide channels. Let 
h{x^ y) be the normalized profile of the waveguide core, 
which is typically assumed to be circularly-shaped with a 
super-Gaussian profile of radius a, (x^, yi) the position of 
the z-th waveguide of the array, and An^ its index change 
from the cladding (substrate) region. One can thus write: 

V{x, y)^-^ Anih{x -Xi,y- yi). (22) 

i 

In the tight-binding approximation, the electric field en- 
velope can be expanded as a superposition of the funda- 
mental guided modes of each waveguide with amplitudes 
that vary slowly along z due to evanescent mode-coupling 
and propagation constant mismatch [29]. In the nearest- 
neighboring approximation, the resulting coupled-mode 
equations, obtained by using e.g. a variational proce- 
dure [30], are precisely of the form (17-19), provided that 
cross-coupHngs among waveguides are negligible and the 
temporal variable t in the quantum HH model is replaced 
by the spatial propagation distance z. If needed, to avoid 
cross-couplings suitable high-index waveguides (or off- 
resonance resonators) could be inserted into the struc- 
ture, as shown in Fig. 1(a). The coupling constants n^^ pn 
entering in the equations (17-19) are basically determined 
by the waveguide separation via an overlapping integral 
of the evanescent fields [sO] , whereas the propagation con- 
stant shifts U + nVt and E± + nVt can be controlled by the 
core index changes An^ of the various waveguides. Let 
us consider, as an example, propagation of light waves 
at wavelength A = 633 nm (red) in waveguides written 
by femtosecond laser pulses in fused silica [16], for which 
Us — 1.45. The normalized core profile of the waveguides 
is assumed to be super-Gaussian-shaped with a core ra- 
dius a = 2.5 /im, namely h{x,y) = exp[— (x^ + y^)^/a^]; 
the index change of the various waveguides is taken to 
be Aui = An + Jn^, where An = 2 x 10~^ is a refer- 
ence index change and Sni a small perturbation, as com- 
pared to An, needed to realize the propagation constant 
shifts. For such waveguides, the coupling constant tz{d) 
(in units of mm~^) between two waveguides placed at 
a distance d (in units of pm) turns out to be very well 
approximated, for distances d >^ 6/im, by the exponen- 
tial relation [3l| k = Aexp(— 7(i), with A 24.6 mm~^ 
and 7 0.466 pm~^. Similarly, the propagation con- 
stant shift a{Sn) of the fundamental waveguide mode 
induced by an index change Jn, added to An, is very 
well approximated (for \Sn\ <C An) by the relation 
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FIG. 1. Schematic of the optical networks that reahze analog 
simulators of the two sets of equations (17-19) for n even 
[Fig. 1(a)] and n odd [Fig. 1(b)]. The filled circles indicate 
either optical waveguides or optical cavities (resonators) that 
are coupled by weak evanescent field. The optical elements 
depicted in Fig. 1(a) by dotted circles indicate high- index core 
waveguides (or off-resonance cavities) that might be inserted, 
if needed, to avoid cross coupling among the optical elements 
in the network. 



(j{5n) —1.7Sn/X. Using such relations, one can thus 
design a waveguide array, with controlled waveguide sep- 
aration and index changes, that realizes the desired cou- 
pling rates a^^, Pn [given by Eq.(20)] and propagation 
constant shifts U + nQ and E± -\- nCl. The total number 
of waveguides needed to simulate the HH model is basi- 
cally determined by the largest phonon number n that 
will be excited in the dynamics. 




propagation distance z (mm) 



FIG. 2. Optical simulations of the HH model in Fock space 
for parameter values ^ = 0.1 mm~^, t = 0.1 mm~^, Q — 
0.1 mm~^, and [7 = (uncorrelated electrons), (a) Refractive 
index profile n{x,y) — ris of the double-chain zig-zag waveg- 
uide array, (b-d) Dynamical evolution of the system corre- 
sponding to the initial state |'0(O)) = \—)®\^)ph- In (b) a 
snapshot of the intensity light pattern \(j){x^y^ z)\^ is shown 
for a few propagation distances 2;, whereas in (c) a snap- 
shot of / dy\(l){x ^ y ^ z)\^ in the (x,^) plane is depicted. Fig- 
ure (d) shows the behavior of the return probability P{z) — 
|('0(O)|'0(2;))|^ = 1/0(2;)!^ numerically obtained from the so- 
lution of the optical Schrodinger equation (solid curve) and 
from the HH model (thin solid curve). 



B. Numerical results 

In this section we provide a few examples of photonic 
simulations of the temporal dynamics of the HH Hamil- 
tonian in Fock space, based on light transport in waveg- 
uide arrays engineered following the procedure outlined 
in the previous subsection, and compare the numerical 
results obtained by the full numerical analysis the op- 
tical Schrodinger equation (21) with those predicted by 
the coupled- mode equations (17-19). As an initial con- 
dition, we will assume the system to be in the state 
— \~) ® corresponding to excitation of zero 

phonon modes and to the two electrons being in the en- 
tangled state |— ) = {l/y/2){c\ ^ — c\ ^4 |)|0), i.e. the 
two electrons occupying the same (either left or right) 
potential well with opposite spins. For such an initial 
state, according to Eq.(16) the initial conditions for the 
Fock-state amplitudes are given by /n(0) = (5^,0 and 
^n(O) = ^n(O) = 0- Such an initial condition simply cor- 
responds to the excitation of the left boundary waveguide 
in the semi- array of Fig. 1(a). The fractional light powers 



trapped in the various waveguides of the array at dif- 
ferent propagation distances z provide a direct mapping 
of the temporal evolution of the occupation probabilities 
l/nP, \9n\^ and \hn\^ of the quantum HH Hamiltonian in 
Fock space. To monitor the temporal dynamics of the 
HH model, we mainly use as a dynamical variable the 
return probability P(t) = |(V^(0)|V^(t))|2 = |/o(t)p. In 
our optical simulator, such a variable simply corresponds 
to the fractional light power that remains in the initially 
excited boundary waveguide. In the simulations, we as- 
sumed constant values of e — ph coupling ^ = 0.1 mm~^ 
and hopping rate t = 0.1 mm~^, whereas a few values 
of the on-site e — e Coulomb energy U are considered, 
corresponding to three different correlation regimes: ab- 
sence of correlation {U = 0), moderate correlation regime 
{U comparable with the hopping rate t); and strong- 
correlation regime {U much larger than t) [32| . 

Figure 2 shows the results obtained for uncorrelated 
electrons {U = 0) and for an oscillator frequency Q = 
0.1 mm~^. In this case, one has cosO = sin^ = 
1/a/2, and thus K,ri = Pn- The refractive index pro- 
file [n(x, y) — Us] of the waveguide array is depicted in 
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FIG. 3. Same as Fig. 2, but in the intermediate correlation 
regime {U = 0.2 mm""*^). 
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FIG. 4. Same as Fig. 2, but in the strong correlation regime 
{U = 10 mm~^). The oscillation frequency is now Q = 
0.2 mm"^ 



Fig. 2 (a). The array comprises 16 waveguides arranged 
in two zig-zag chains with non-uniform spacings. The 
evolution of the light intensity along the array is shown 
in Figs. 2(b) and (c). In particular, in Fig. 2(b) a snap- 
shot of the intensity light distribution |^(x,?/,z)p is de- 
picted for a few propagation distances (up to z = 8 
cm) as obtained by numerical simulations of the opti- 
cal Schrodinger equation (21), whereas in Fig. 2(c) the 
snapshot of J dy\(j){x,y^ z)]"^ in the (x^z) plane is shown. 
Note that a relatively small number of waveguides (~ 13) 
are excited in the dynamics. The corresponding behav- 
ior of the return probability P{z) is depicted in Fig. 2(d) 
(solid curve) and compared with the behavior obtained 
by direct numerical analysis of the HH equations (17-19) 
(thin solid curve). Note that, even though in our optical 
network high-index waveguides to avoid weak cross cou- 
plings have not been inserted, the agreement with the 
HH model [Eqs. (17-19)] is rather satisfactory. Figure 3 
shows similar results obtained in the moderate correlated 
regime {U = 0.2 mm~^), corresponding to cos^ 0.851, 
sin^ 0.526. Note that in this case hCn Pm and thus 
the two zig-zag chains forming the array are noticeably 
asymmetric. 

The strong correlation regime is attained when the on- 
site interaction strength U is much larger than the hop- 
ping amplitude t, corresponding to cos^ 1 and sin^ :^ 
0. In this case ~ and thus the phonon modes are 
coupled to the electronic degree of freedoms solely via 
the two dressed states |— ) and — 1+), whose en- 

ergy levels are nearly degenerate and equal to U and 
Ej^ '^U{1 +4t^//7^), respectively. In Fock space, the 



e — ph couphng dynamics is thus reduced to the cou- 
pled equations for the occupation amplitudes fn and 
namely [see Eqs. (17) and (18)] 

= {U ^ Vtn)fn - I^n9n-l " I^n+l9n+l (23) 
= (^+ + n^)9n - l^nfn-1 - ^n+l/n+1 (24) 

where kn = V2ng. Correspondingly, the optical analog 
simulator reduces to a simple linear chain of waveguides 
with non-uniform spacing, as shown in Fig. 4(a). An ex- 
ample of optical simulation of the HH model in the strong 
correlation regime is shown in Fig. 4 for parameter values 
g = 0.1 mm~^, t = 0.1 mm~^, Q = 0.2 mm~^ and 
U = 10 mm~^. Note that the dynamics of occupation 
amplitudes shows in this case a typical periodic behav- 
ior with a period given by 27r/r^. In our optical simula- 
tor, the periodic revival is clearly visualized by the self- 
imaging property of the waveguide array, in which the 
initial light periodically returns into the initially-excited 
waveguide [see Fig. 4(c)]. Such a periodic behavior can be 
explained after observing that, in the strong correlation 
regime ^ oo, the phonon modes \n)ph couple solely 
to two dressed electronic states |— ) and |+). Neglecting 
the small energy shift {E-^ — U) of dressed levels |— ) and 
|+), the energy levels of the e — ph coupled system [the 
eigenvalues of Eqs. (23) and (24) with = [/] is analo- 
gous to a Wannier- Stark ladder with energy spacing ft, 
namely one has Ei = U ^ IQ - 2g^ jVi (/ = 0, 1, 2, 3, ...). 
In fact, in this limit the HH model is equivalent to a 
two-site tunneling system strongly coupled to a quantum 
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harmonic oscillator, such as the two-site Holstein model 
in atomic (strong e — ph coupling) limit [33] , describ- 
ing the hopping dynamics of one electron in a two-site 
potential interacting with local phonon modes, or the 
Jaynes-Cummings model of quantum optics, describing 
the interaction of a two-level atom with a quantized os- 
cillator in the so-called deep strong coupling regime [34] • 
In the atomic (strong-coupling) limit, i.e. for a vanish- 
ing hopping rate, the Holstein model is exactly integrable 
using a Lang-Firsov transformation 0, [35| , yielding an 
equally-spaced discrete energy spectrum corresponding 
to small Holsten polaron modes . In our optical realiza- 
tion, the equally-spaced energy spectrum is responsible 
for the periodic self-imaging of the optical beam as it 
propagates along the array, with a spatial period 27r/Q. 
In our realization of the HH model in Fock space, the 
periodic self-imaging effect can be explained in terms of 
generalized Block oscillations of a single particle in a non- 
homogeneous tight-binding lattice with a dc bias. In fact, 
the tight-binding lattice realization of the HH model in 
the strong correlation regime, as described by Eqs.(23) 
and (24) with E-^ /7, belongs to a class of exactly- 
solvable lattice models which shows an equally-spaced 
Wannier-Stark ladder spectrum, as shown in Ref.[36,]. In 
this lattice model, the dc biased in represented by the lin- 
ear gradient terms nQ in Eqs.(23) and (24), and Bloch os- 
cillations are not washed out in spite of lattice truncation 
and non- homogeneity of hopping rates Ref . (36^ . . If the 
small energy shift between the dressed energy states |— ) 
and 1+) is taken into account, the revival (self-imaging) 
dynamics would be only approximate because the energy 
levels of H are no more exactly equally spaced (see, for 
instance, Q). However, for the parameter values used 
in the simulation of Fig. 4, washing out of the self-imaging 
would be visible after several self-imaging periods, and it 
is thus not observed in the figure. 

IV. CONCLUSIONS 

In this work, a classical analog simulator of the 
Hubbard-Holstein model, describing electron-phonon 
coupling in presence of electron correlations, has been 
theoretically proposed based on light transport in engi- 
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